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Abstract. Binary black hole simulations starting from quasi-circular (i.e., zero radial 
velocity) initial data have orbits with small but non-zero orbital eccentricities. In 
this paper the quasi-equilibrium initial-data method is extended to allow non-zero 
radial velocities to be specified in binary black hole initial data. New low-eccentricity 
initial data are obtained by adjusting the orbital frequency and radial velocities to 
minimize the orbital eccentricity, and the resulting (~ 5 orbit) evolutions are compared 
with those of quasi-circular initial data. Evolutions of the quasi-circular data clearly 
show eccentric orbits, with eccentricity that decays over time. The precise decay 
rate depends on the definition of eccentricity; if defined in terms of variations in the 
orbital frequency, the decay rate agrees well with the prediction of Peters (1964). The 
gravitational waveforms, which contain ~ 8 cycles in the dominant I — m — 2 mode, 
are largely unaffected by the eccentricity of the quasi-circular initial data. The overlap 
between the dominant mode in the quasi-circular evolution and the same mode in the 
low-eccentricity evolution is about 0.99. 
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1. Introduction 

The inspiral and merger of binary black holes is one of the most promising sources for 
current and future generations of interferometric gravitational wave detectors such as 
LIGO and VIRGO [1, 2]. The initial LIGO detectors, which are currently operating 
at design sensitivity, could detect binary black hole inspirals up to distances of several 
hundred megaparsecs. In order to take full advantage of the sensitivity of these detectors, 
detailed knowledge of the gravitational waveform is required. 

Recent breakthroughs in numerical relativity have allowed several research groups 
to simulate binary black hole inspirals for multiple orbits [3, 4, 5, 6, 7]. Because of 
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Figure 1. Evolution of quasi-circular initial data. The left panel shows the proper 
separation s between the apparent horizons, computed at constant coordinate time 
along the coordinate line connecting the centers of the horizons, and the right panel 
shows its time derivative ds/dt. This evolution was run at three different resolutions, 
with the medium and high resolution tracking each other very closely through the run. 



the large computational cost of these simulations, only a small number of orbits can 
be followed. Therefore it is important to begin these simulations with initial data that 
closely approximate a snapshot of a binary black hole system that is only a few orbits 
from merger. During the inspiral, the orbits of binary compact objects circularize via 
the emission of gravitational waves [8] , so binaries formed from stellar evolution (rather 
than dynamical capture) are expected to have very small eccentricities by the time they 
enter the sensitive band of ground based detectors. Because of this, the assumption of a 
quasi-circular orbit (i.e., zero radial velocity) has been widely used in the construction 
of binary black hole initial data [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. 
Specifically, quasi-equilibrium data [17] and the "QC-sequence" [24] of puncture 
data [25] seem to be the most popular, and both of these assume a quasi-circular orbit. 
However, inspiraling compact objects have a small inward radial velocity, and neglecting 
this velocity when constructing initial data will lead to eccentricity in the subsequent 
evolution, as discussed in the context of post-Newtonian theory in Ref. [26], and found 
numerically in Refs. [27]. 

The Caltech/Cornell collaboration has recently completed successful long-term 
simulations of inspiraling binary black holes [6] using a pseudo-spectral multi-domain 
method. This technique was used to evolve a particular quasi-circular quasi-equilibrium 
binary black hole initial data set (coordinate separation d = 20 from Table IV of 
Ref. [17]). Figure 1 shows the proper separation s between the horizons and the 
radial velocity ds/dt as functions of time for this evolution. The rapid convergence 
afforded by spectral methods is apparent; the medium and high resolutions are nearly 
indistinguishable on the plot. Eccentricity of the orbit in the form of oscillatory 
variations in s and ds/dt is, unfortunately, also clearly apparent. 
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This noticeable eccentricity suggests two questions: First, how can initial data with 
the appropriate black hole radial velocities be constructed for non-eccentric inspirals? 
Second, how do evolutions of quasi-circular initial data differ from those of non-eccentric 
initial data? This paper addresses both questions. In Sec. 2, we incorporate nonzero 
radial velocities into the quasi-equilibrium method to construct binary black hole initial 
data. This results in one additional parameter for equal mass initial data, the radial 
velocity v r . Section 3 briefly discusses our numerical methods. Section 4 describes 
how we choose v r and the the orbital frequency VLq for equal mass co-rotating binary 
black holes, and presents numerical evolutions of the resulting low-eccentricity initial 
data. This section also presents convergence tests of these binary black hole evolutions; 
we examine both convergence with respect to spatial resolution and convergence with 
respect to the radius of the outer boundary of the computational domain. Section 5 
examines the differences between evolutions of quasi- circular initial data and low- 
eccentricity initial data. We close with a summary and discussion of these results in 
Sec. 6. 

2. Quasi-equilibrium data with nonzero radial velocity 

In this section we extend the quasi-equilibrium approach [14, 16, 17, 20] to allow 
specification of nonzero radial velocities of the black holes. We proceed in three 
steps: First, we summarize the construction of quasi-equilibrium data using co-rotating 
coordinates [17, 20]. Second, we show that the identical quasi-circular initial data can 
be obtained by solving essentially the same equations but in an asymptotically inertial 
coordinate system; the major difference is that one must require the black holes to move 
on circular trajectories, rather than remaining fixed in the coordinate system. Third, 
we generalize from black holes moving on circular trajectories to black holes moving on 
inspiral trajectories. 

2.1. Overview 

We use the nomenclature of Ref. [17]; the spacetime line element is written in the usual 
3 + 1-form, 



where 7^ is the 3-metric induced on a t = constant spatial hypersurface, a is the 
lapse function and f3 l is the shift vector. Latin indices label spatial coordinates, and 
Greek indices label spacetime coordinates. The extrinsic curvature of the hypersurface 
is defined by 



to the sliced . We use the extended conformal thin sandwich formalism [28, 29] to 

\ Since is a spatial tensor, K^ v n v = 0, its spatial components Kij carry all its information. Almost 
all tensors in this paper are spatial, and we use spatial indices here whenever possible. 



ds 2 = -a 2 dt 2 + ^(fa? + p i dt)(dx j + (3 j dt), 



(1) 
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construct constraint-satisfying initial data. In this approach, the three-dimensional 
metric is spit into a conformal metric %j and a positive conformal factor ip, 

lij = (3) 
and the extrinsic curvature is split into trace and trace-free parts, 

Kn = Aij + --YijK. (4) 

The freely specifiable data consist of the conformal metric 7^, its time derivative 
Uij = dt^ij (which is taken to be trace free), the mean curvature K = Kiff 3 , and 
its time derivative d t K. It follows that the trace-free part of the extrinsic curvature 
takes the form 

An = — [{Lp)ij - V>%] = ij- 2 A tv Aij = [(£/%• - u tj \ , (5) 



where 



(L/3)« = 2V^> _ H 7 «v fc /3 fc , (lp) ij = 2V^'> - ^V k (3 k . (6) 



The symbols (L/3) u and (L/?)*- 7 represent the conformal Killing operators in physical 
and conformal space, respectively, and are related by (L/3) 1 - 7 = ■0 _4 (L/9) y . Indices 
on conformal tensors are raised and lowered with the conformal metric, for example, 
— Iiklji(] u l3) kl — ip 4 (L/3)jj. Furthermore, V and V denote the physical and 
conformal spatial covariant derivative operators, and the conformal lapse is defined by 
a = tp 6 a. Inverting Eq. (5) yields 

= drfij = -2aAij + (L/3)ij. (7) 

Substituting these relations into the constraint equations and into the evolution 
equation for the extrinsic curvature, one arrives at a system of five elliptic equations, 
often referred to as the extended conformal thin sandwich (XCTS) equations: 

V 2 ^ - \W - -^K 2 ^ + ^A^Aij = 0, (8a) 
V, - - % (^fi«) = 0, (86) 

V 2 (a^ 7 ) - (a^j 7 ) [f + l-^+U-^Ai ,] = ~^{d t K - (3 k d k K). (8c) 

L a iZ 8 J 

Here R denotes the trace of the Ricci tensor of 7^. These equations are to be solved 
for ip, a and f3 l ; given a solution, the physical initial data ( , jij,Kij) are obtained from 
Eqs. (3)-(5). 

Note that a solution of the XCTS equations includes a shift vector f3 l and a lapse 
function a = ip 6 a. If these values of lapse and shift are used in an evolution of the 
constructed initial data, then the time derivative of the mean curvature will initially 
equal the freely specifiable quantity d t K, and the trace-free part of the time derivative 
of the metric will initially equal ip 4 Uij. Thus, the free data of the XCTS equations allow 
direct control of certain time derivatives in the evolution of the initial data. 
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The next step is to choose the free data that correspond to the desired physical 
configuration. The quasi-equilibrium quasi-circular orbit method of constructing binary 
black holes [17, 20] (see also Refs. [11, 12, 14]) provides a framework for many of these 
choices. This method is based on the fact that the inspiral time scale for a binary 
compact object is much larger than the orbital time scale, so that time derivatives 
should be very small in the co-rotating coordinate system. Furthermore, the black holes 
should be in equilibrium, which provides conditions on the expansion 9 and shear of 
the outgoing null geodesies passing through the horizon. The complete set of physically 
motivated choices for the free data within the quasi-equilibrium method are 

Uij = 0, (9 a) 

d t K = 0, (96) 

tp — > 1, a — > 1, as r — > oo, (9c) 

ft -> (n x r)*, as r -> oo, (9(f) 

d t is tangent to <Sah, (9 e ) 

6 = on <S, (9f) 

Uij = on S, (9g) 

where S denotes the location of the apparent horizons in the initial data surface, and 
Sah is the world tube of the apparent horizon obtained by evolving the initial data 
with lapse a and shift ft. The first two conditions are the assumptions that the time 
derivatives are small. The boundary conditions in Eqs. (9 c) and (9(f) enforce asymptotic 
flatness and co- rotation. The orbital frequency f2 entering Eq. (9(f) can be chosen by the 
effective potential method [9] or the Komar-mass ansatz [11], with similar results [20]. 

To discuss the remaining conditions, we need to introduce a few additional 
geometrical quantities. Denote by s l and s l the physical and conformal outward-pointing 
spatial unit normals to S. They obey the relations 

sV 7o - = l, sV'7tf = l, s*^" 2 ?. (10) 

Then introduce the induced metric on S in physical and conformal space by hij = 
1ij — SiSj, and = respectively. Because n^s^O, the space-time components 

of the unit normal are given by s M = [0, s 1 ]. The outward-pointing null normal to S can 
then be written as 

F = (n" + *") . (11) 

Equation (9e) simply means that the apparent horizon is initially at rest when the 
initial data is evolved in the co-rotating coordinate system. It implies that the shift 
must take the form 

ft = aS i + ft on S, (12) 

where ft is tangent to S. Equation (9f) ensures that S is an apparent horizon, and 
implies a boundary condition on the conformal factor, 



8a 



ip 7ii ~„ 1 



~s k d k ^j = § 1 <P (hp^j-Uij -^-h i:> V i s j + -Ki , 6 . (13) 



4 —J 6 - 
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Finally, Eq. (9g) — which forces the apparent horizon to be in equilibrium — restricts (3^ 
to be a conformal Killing vector within the surface S, 

(IsW = 2L>^ ) - ~tiW k p* = 0, (14) 

where Di is the covariant derivative compatible with h^. As discussed in detail in 
Refs. [17, 20], controls the spin of the black holes in addition to the spin required for 
co-rotation. 

Quasi-equilibrium considerations have now led us to choices for half of the free 
data (Hij and d t K) for the XCTS equations, and for all boundary conditions except a 
lapse boundary condition on the horizon S. As argued in Ref. [17], Eqs. (9a)-(9e) are 
compatible with any spin of the black holes, with any choice of boundary conditions for 
the lapse on S, and with any choice of ^ and K. For concreteness, we choose 

la = fa, ( 15a ) 
K = 0, (15b) 

d r (at/j) = on S, (15c) 

where is the Euclidean metric. The last two conditions, Eqs. (156) and (15c), are 
gauge choices [17]. The choice of the conformal metric, however, does influence the 
physical gravitational radiation degrees of freedom of the system. Since a black hole 
binary is not conformally flat at second post-Newtonian order [30], our simple choice of 
conformal flatness, Eq. (15a), is probably responsible for the initial burst of unphysical 
gravitational radiation found in the evolution of these initial data. 



2.2. Initial data in an asymptotically inertial frame 

It is possible to re-formulate the quasi-equilibrium method in asymptotically inertial 
coordinates in such a way that identical physical initial data are obtained. To do so, we 
solve the XCTS Eqs. (8a)-(8c) with the same choices for the free data and boundary 
conditions, except that Eqs. (9d) and (9e) are replaced by 

f3 l — > as r — > oo, (16a) 

d t + CoA is tangent to S AH , where Cot = (^o x r) 1 . (166) 

The second condition implies that the apparent horizons move initially with velocity 
£r Qt , i.e. tangent to circular orbit trajectories. 

Let (ip co , P l co ,a co ) be the solution to the XCTS equations in the co-rotating 
coordinates. We show in Appendix A that the solution in the asymptotically inertial 
coordinates is (ip, a) = (ip co , Pl Q — Cot> a co), and that this solution leads to the same 
physical metric 7^ and extrinsic curvature Kij as the original solution in co-rotating 
coordinates. The proof of this relies on two observations: First, the shift enters the 
XCTS equations and the boundary conditions (almost) solely through the conformal 
Killing operator, (L/3) 3 ; and second, Q ot is a conformal Killing vector, (L^ot) 1 - 7 = 0, for 
the conformally flat case considered here. Hence the term — ^* ot that is added to f3 l co 
drops out of the equations. 
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In Appendix A, we also show that Eq. (166) and the shear condition Eq. (9g) require 
the shift on the inner boundary S to take the form 

(? = as i -g ot + ( i on S, (17) 

where ( l is a vector that must be tangent to S (C s i — 0) an d must be a conformal 
Killing vector within the surface S: 

a lj = & = (l s Cy'. (18) 

Comparing Eq. (17) with Eq. (12), we see that the vector ( l plays the role of /3l in 
the earlier treatment; choosing it as a rotation within S will impart additional spin to 
the black holes in addition to co-rotation, as described in detail in Ref. [20]. Note that 
at large radii the comoving shift (3 l co is a pure rotation, since (3 l co — >■ £* ot [Eq. (16a)] and 
V J £rot 1S antisymmetric [Eq. (166)]. 

2.3. Initial data with nonzero radial velocity 

After rewriting the standard quasi-equilibrium method in an asymptotically inertial 
frame, it is straightforward to incorporate nonzero initial radial velocities for the black 
holes. As discussed in Sec. 2.2, quasi-circular initial data can be generated by specifying 
that the horizons move initially on circles in an asymptotically inertial coordinate 
system. This is accomplished by the shift boundary conditions in Eqs. (16a) and (166). 
We include initial radial velocities simply by requiring the black holes to move initially 
on inspiral rather than circular trajectories. 

Consider the problem of giving a black hole located a distance r from the origin 
an initial radial velocity v r . This can easily be accomplished by replacing the boundary 
conditions in Eqs. (16a) and (166) with 

ft -> as r -> oo, (19a) 

9t + £ins P <% is tangent to S AU , where £- nsp = (fl x r) 4 + vj—. (196) 

As before, we place the center of rotation at the origin of the coordinate system. Note 
that £? nsp is still a conformal Killing vector, (L^p) 1 - 7 = 0, for the conformally flat case 
considered here. Therefore the analysis in Appendix A of the boundary conditions in 
Eqs. (166) and (9a) also applies to Eqs. (196) and (9a), and so we find that the inner 
shift boundary condition must be of the form 

ft = as' - £Lp + C, on 5, (20) 

where ( l is a conformal Killing vector within S. 

The boundary conditions in Eqs. (19a) and (196) depend on two parameters, the 
orbital frequency Qo and a radial velocity v r (or, more precisely, an overall expansion 
factor v r /ro, reminiscent of the Hubble constant). For unequal mass binary systems 
the needed radial velocities for each hole would be different, but the needed expansion 
factors, v r /r Q , are expected to be the same for the two holes. 
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The changes discussed in Sec. 2.2 are superficially similar to the changes discussed 
in Sec. 2.3, yet the former amounts to a mere coordinate transformation while the 
latter produces different physical initial data. This can be understood by noting that 
the change from co-rotating coordinates [Eqs. (9d) and (9e)] to inertial coordinates 
[Eqs. (16a) and (166)] is accomplished by adding the same conformal Killing vector field 
£* ot to the shift at both inner and outer boundaries, but the change from Eqs. (9(f) and 
(9e) to initial data with nonzero radial velocity [Eqs. (19a) and (196)] is accomplished 
by adding different conformal Killing fields to the shift on different boundaries: Q ot at 
the outer boundary and Q nsp at the inner boundaries. Only in the former case can the 
change be expressed as a global transformation of the shift of the form f3 % — > (3 l + Q ot . 

3. Numerical methods 

The initial value equations are solved with the pseudo-spectral elliptic solver described in 
Ref. [31]. This elliptic solver has been updated to share the more advanced infrastructure 
of our evolution code and is now capable of handling cylindrical subdomains. This 
increases its efficiency by about a factor of three over the results described in Ref. [31] 
for binary black hole initial data. 

The Einstein evolution equations are solved with the pseudo-spectral evolution code 
described in Ref. [6]. This code evolves a first-order representation [32] of the generalized 
harmonic system [33, 34]. We use boundary conditions [32] designed to prevent the influx 
of unphysical constraint violations and undesired incoming gravitational radiation, while 
allowing the outgoing gravitational radiation to pass freely through the boundary. The 
code uses a fairly complicated domain decomposition. Each black hole is surrounded by 
three concentric spherical shells, with the inner boundary of the inner shell just inside 
the horizon. The inner shells overlap a structure of 24 touching cylinders, which in turn 
overlap a set of outer spherical shells — centered at the origin — which extend to large 
outer radius. Outer boundary conditions are imposed only on the outer surface of the 
largest outer spherical shell. We vary the location of the outer boundary by adding more 
shells at the outer edge. Since all outer shells have the same angular resolution, the cost 
of placing the outer boundary farther away (at full resolution) increases only linearly 
with the radius of the boundary. Some of the details of the domain decompositions used 
for the simulations presented here are given in Table 1. 

4. Choice of orbital frequency and radial velocity 

We now describe how to construct binary black hole initial data sets with low orbital 
eccentricity. This is done by tuning the freely adjustable orbital parameters f2 and v r 
iteratively to reduce the eccentricity of the inspiral trajectories. For each iteration we 
choose trial orbital parameters Q Q and v r , evolve the corresponding initial data, analyze 
the resulting trajectories of the black holes, and update the orbital parameters to reduce 
any oscillatory behavior in quantities like the coordinate separation of the black holes 
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d(t), the proper separation between the horizons s(t), or the orbital frequency cu(t). All 
of these quantities (and many others) exhibit similar oscillatory behavior; we choose d(t) 
as our primary diagnostic during the tuning process because it is most easily accessible 
during the evolutions. 

To make this procedure quite explicit, we begin by evolving quasi-circular initial 
data for about two orbits. Then we measure the time derivative of the measured 
coordinate separation of the holes d(t) (in the asymptotic inertial coordinates used 
in our code [6]) as illustrated for example in Fig. 3. We fit this measured d(t) to a 
function of the form: 

d{t) = A + Att + Bsin(Lut + if), (21) 

where A , A 1: B, u, and (p are constants determined by the fit. The A + A x t part 
of the solution represents the smooth inspiral, while the Bsm(cut + ip) part represents 
the unwanted oscillations due to the eccentricity of the orbit. For a nearly circular 
Newtonian orbit, B is related to the eccentricity e of the orbit by e = B/uod. So 
reducing the orbital eccentricity is equivalent to reducing B. The values of the orbital 
parameters Qq and v r are now adjusted iteratively to make the coefficient B in this fit 
as small as desired. After each adjustment of Qq and v r , the initial value equations 
described in Sec. 2 [in particular, using the boundary condition (196) which depends on 
Qq and v r ] are solved completely (to the level of numerical truncation error). 

For this paper, our goal is to reduce B, and hence the orbital eccentricity, by about a 
factor of ten compared to quasi-circular initial data. This level of reduction is sufficient 
to allow us to evaluate the significance of the orbital eccentricity inherent in quasi- 
circular initial data. A variety of methods could be used to find orbital paramters that 
make B small. One possibility is simply to evaluate B(Qo,v r ) numerically as described 
above, and then to use standard numerical methods to solve the equation B(Q , v r ) = 0. 
Since our goal in this paper is to reduce B by about a factor of ten, simple bisection 
root finding methods are sufficient. 

A more efficient method is to use our knowledge of the behavior of nearly circular 
orbits to make informed estimates of the needed adjustments in the orbital parameters. 
Evaluating the fit Eq. (21) at the initial time t — 0, we see that the ellipticity-related 
component B sin(u;t + <p) contributes _Bsin(<^)/2 to the radial velocity of each hole and 
Bio cos(<^)/2 to its radial acceleration. (The factor 1/2 arises because d measures the 
distance between the holes.) For a Newtonian binary, this eccentricity- induced radial 
velocity can be completely removed by changing the initial radial velocity by 

SVr = ( 22) 

Furthermore, changing the orbital frequency Qo by a small amount 5Qo changes the 
radial acceleration of each black hole by the amount QoSQodo, where do = d(0) is the 
initial separation of the holes. Thus the change <5f2 needed to remove the eccentricity- 
induced initial radial acceleration, Bu cos(f)/2, is 

SQ = BuJCOS (<P) w Bcos(ip) ^ 
2c?of2o 2g?o 
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Table 1. Summary of evolutions presented in this paper. The labels 'QC', 'E' and 'F' 
refer to the different initial data sets, with numerical suffix ('El', 'E2' etc.) denoting 
different values of the initial outer boundary radius of the evolutions, i? ou ter- 



Label 


Initial data 


m a "dm ^ outer a PP r ox # of points 
shells low med. high 


QC 


M A dm^o = 0.029792, v r = 0.0 
Jadm/M adm = 0.98549 
Mi„/M AD M = 0.50535 


133 8 52 3 64 3 76 3 


El 
E2 


Madm^o = 0.029961, v r = -0.0017 
Jadm/M adm = 0.99172 
M irr /M AD M = 0.50524 


171 10 59 3 66 3 74 3 
293 18 64 3 72 3 81 3 


Fl 
F2 
F3 


M AD mO = 0.029963, v r = -0.0015 
Jadm/M adm = 0.99164 
M irr /M AD M = 0.50525 


133 8 52 3 64 3 76 3 
190 12 55 3 66 3 78 3 
419 28 62 3 74 3 87 3 



Equations (22) and (23) still hold approximately for relativistic binaries. We have found 
that simultanously adjusting v r and Qq by Eqs. (22) and (23) typically reduces B by 
about a factor of ten. 

The smallest eccentricity data set produced here (by the simple bisection method 
described above) is labeled 'F', and the data from the next to last iteration of this 
method is labeled 'E'. These initial data sets, together with the quasi-circular data 
labeled 'QC were evolved with multiple numerical resolutions and with multiple outer 
boundary locations; Table 1 summarizes these evolutions. The orbital frequency used 
in the final evolution is only 0.6 per cent larger than the value of VL used in the quasi- 
circular case. As expected, this change is comparible to the magnitude of the radial 
velocity v r in the low eccentricity case. The smallness of these quantities shows that the 
quasi-circular approximation is quite good. 

Figure 2 shows the orbital phase (as measured by the coordinate locations of the 
centers of the apparent horizons) for the evolutions of quasi-circular initial data, QC, 
and the least-eccentric initial data, Fl, F2, and F3. (The numerical suffix, Fl, F2, etc., 
denotes simulations with different values of the outer boundary radius as defined in 
Table 1.) These evolutions proceed for about five orbits and then crash shortly before 
the black holes merge. The upper left inset shows differences between the orbital phase 
computed with different resolutions for the QC and the F2 runs. The phase difference 
between the high and low resolution runs is < 0.35 radians, which is a good estimate of 
the error in the low resolution run. The phase difference between the medium and high 
resolution runs drops to ~ 0.02 radians, which can be taken as the error in the medium 
resolution run. Between low and medium resolutions, the error drops by about a factor 
of 20. Assuming exponential convergence, the error of the high resolution run should 
be smaller by yet another factor of ~ 20, i.e. < 0.001 radians. The lower right inset 
in Fig. 2 shows phase differences between evolutions of the same initial data, but run 
with different outer boundary radii. These differences are small, so we do not expect 
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Figure 2. Evolution of the orbital phase. The main panel shows the phase of the 
trajectories of the centers of the apparent horizons as a function of time for the quasi- 
circular (dotted curves) and low-eccentricity (solid curves) initial data. The top left 
inset shows the phase differences between different resolution runs, which decreases 
at higher resolutions. The lower right inset shows the difference in the orbital phase 
between evolutions with different outer boundary locations. 
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Figure 3. Radial velocity during evolutions of quasi-circular and low-eccentricity 
initial data. The left panel shows the coordinate velocity d(t), the right panel the 
velocity determined from the intra-horizon proper separation s(t). 



the influence of the outer boundary on our results to be significant. Our analysis in 
Sec. 5 is based mostly on comparisons between the high resolution QC and F2 runs. 

Figure 3 illustrates the radial velocities (determined from the time derivatives of 
both the coordinate and the intra-horizon proper separations) for the quasi-circular run 
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QC and for the two low-eccentricity runs E and F. Orbital eccentricity causes periodic 
oscillations in these curves; the amplitudes of these oscillations are clearly much smaller 
in runs E and F than in run QC. By fitting the proper separation speed ds/dt to a linear 
function plus sinusoid, ds/dt = Aq + A\t + B sin(W + tp), the approximate amplitude 
of the oscillations can be estimated. We find -Bqc ~ 5.5 • 1CT 3 , Be ~ 5.8 • 10~ 4 , and 
Bp ps 4.1 • 10~ 4 . This confirms that we have succeeded in our goal of reducing the 
oscillations by an order of magnitude. These fits are not very accurate because the fit 
must cover at least one period of the oscillations, and significant orbital evolution occurs 
during this time. If we vary the fit interval 40 < t/M ADM < T by choosing T between 
300 and 450, the quoted amplitudes Aq C ,e,f change at about the 10% level. 

The coordinate separation d(d)/dt shows some noise at early times as the binary 
system equilibrates and an initial burst of 'junk' gravitational radiation travels outward. 
There are also short-lived, high-frequency features apparent in Fig. 3 at intermediate 
times. The earlier feature occurs at £/Madm ~ 140 for the QC run, £/Madm ~ 200 for 
F2, and t/MADM ~ 300 for E2; these times coincide with the light-crossing time to the 
outer boundary. We believe that this early feature is caused by a small mismatch 
between the initial data and the outer boundary conditions used by the evolution 
code; this mismatch produces a pulse that propagates inward from the outer boundary 
starting at t — 0. A later (and larger) feature occurs at t/M ADM ~ 280 for the QC 
run, t /Af adm ~ 400 for F2, and at t/M ADM ~ 600 (off the scale of Fig. 3) for E2. 
This later feature occurs at twice the light-crossing time, and is caused by reflection 
of the initial 'junk' gravitational radiation burst off of the outer boundary. The outer 
boundary conditions used in this paper perform well for the physical gravitational- wave 
degrees of freedom [32], but comparatively poorly for the gauge degrees of freedom (as 
demonstrated in recent tests [35]). These results plus the observation that the high- 
frequency features in Fig. 3 are greatly diminished in less gauge- dependent quantities 
like ds/dt suggests that these features may be caused by perturbations in the gauge or 
coordinate degrees of freedom of the system. 

Figure 4 shows the orbital trajectories of the centers of the black holes during 
evolutions of the low-eccentricity initial data E §, and the quasi-circular initial data QC. 
The low-eccentricity run forms a smooth spiral with no apparent distortion. In contrast, 
the evolution starting from quasi-circular initial data has clearly visible irregularities. 

5. Comparing quasi-circular and low-eccentricity initial data 

Figures 3 and 4 show clearly that evolutions of the quasi-circular initial data, QC, are not 
the same as those of the low-eccentricity initial data, F. In this section, we characterize 
and quantify these differences in more detail. 

§ We plot the evolution El because it was pushed somewhat closer to merger than the F runs; the 
trajectories of the E runs are indistinguishable from those of the F runs on the scale of this figure. 




Figure 4. Trajectories of the center of the apparent horizons in asymptotically inertial 
coordinates for the runs El (left plot) and QC (right plot). The solid/dashed line 
distinguish the two black holes; the circles and ellipsoids in the left figure denote the 
location of the apparent horizon at the beginning and end of the evolution. 



5.1. Time shift 

The black holes approach each other more quickly in the QC run, with the time of 
coalescence appearing to be about 60Madm earlier than in the F2 run. Figure 2, for 
example, shows that the orbital phase increases more quickly during the QC run, with a 
late time phase difference of about tt (almost a full gravitational wave cycle) compared 
to the F2 run. Similar differences are also seen in the graphs of the proper separation 
and orbital frequency shown in the upper panels of Fig. 5. 

We find that most of the difference between the QC and F runs is just a simple 
coordinate time shift. The dashed lines in the upper panels of Fig. 5 represent the QC 
evolution shifted by AT = 59Madm- With this time shift, the QC evolution oscillates 
around the low-eccentricity F2 run. Therefore, the apparent earlier merger time of the 
QC run is just a consequence of the fact that coordinate time t = in the QC run 
represents a later stage in the inspiral than it does in the F2 evolution. The QC and F2 
runs were started with the same spatial coordinate separation at t — 0; however, this 
point is the apocenter of the slightly eccentric QC orbit, so the point in the F2 run with 
the same phase (measured from merger) has smaller separation. 

The lower left panel of Fig. 5 shows the proper separation difference, 5s = 
sp(t) — SQc(t — AT), which emphasizes the oscillations of the QC evolution around 
the F2 orbit. These differences are plotted for three different time shifts AT. The right 
panels of Fig. 5 present information about the orbital angular frequency uj as determined 
from the coordinate locations of the centers of the apparent horizons. The upper right 
panel shows ui for evolutions of QC and F2 initial data. Time-shifting the QC run by 
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Figure 5. Proper separation (left) and orbital frequency (right) for evolutions of the 
QC and F initial data. The lower panels show the differences between the time-shifted 
QC and the F2 runs. The dotted lines in the lower panels show the differences between 
the El and F2 runs, providing an estimate of the remaining eccentricity in the F2 run. 



the same AT = 59Madm also lines up the frequency curves very well. The lower right 
plot shows the difference in orbital frequency between the F2 run and the time-shifted 
QC run, 5u = up(t) — ujQc(t — AT). The differences 5s and 5u are very sensitive to 
the time offset AT applied to the QC run. In particular, at late times, when s and u 
vary rapidly, even a small change in AT causes the differences to deviate significantly 
from their expected oscillatory behavior around zero. Looking at both 5s and 8w, we 
estimate a time offset AT/Madm = 59 ± 1 between the QC run and the F runs. 



5.2. Measuring eccentricity 

The evolution of the F initial data appears to have very low orbital eccentricity, so it can 
be used as a reference from which the eccentricity of the QC run can be estimated. We 
can define an eccentricity for the QC evolution, for example, from the relative proper 
separation, 

e s = ^1, (24) 
s 

where this equation is to be evaluated at the extrema of 5s. Similarly, we can define a 
different measure of eccentricity from the variations in u> or ut by evaluating 
\5u\ 

e " = ^7 

at the extrema of 5u. The factor of two in the definition of arises from angular 
momentum conservation, which makes the orbital frequency proportional to the square 
of the radius of the orbit. In Newtonian gravity, e s = to first order in eccentricity. 
Since the F initial data results in a factor of ten smaller oscillations in ds/dt than the QC 



(25) 
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Figure 6. Orbital eccentricity of the QC evolution estimated from variations in proper 
separation, e s , and from variations in orbital frequency, e w . Also shown in this log- log 
plot are best-fit power laws to each set of data, as well as the scaling predicted by 
Peters [8] with power 19/12 w 1.58. 



data, we expect these eccentricity estimates to be affected by the residual eccentricity 
of the F run at only the 10% level. 

The orbital eccentricity of the QC run, estimated using Eqs. (24) and (25), is 
plotted as a function of proper separation between the black holes in Fig. 6. We see 
that these eccentricities decay during the inspiral, as expected. Within our estimated 
10% errors, these eccentricities are consistent with a power law dependence on the proper 
separation, e oc s p . The eccentricity e s based on the proper separation is consistently 
somewhat larger than e w , and it decays somewhat more slowly. Peters [8] derived the 
evolution of the orbital eccentricity during an inspiral due to the emission of gravitational 
waves using the quadrupole approximation. His result in the e< 1 limit predicts that 
e oc a 19 / 12 , where a is the semi-major axis of the orbit and where the constant of 
proportionality depends on the initial conditions. Using a « s/2, his formula predicts 
that the eccentricity should decay as 

e oc s 19 ' 12 . (26) 

Figure 6 confirms that e w follows this prediction quite closely, while e s has a somewhat 
smaller power law exponent. 

The eccentricities measured here are actually the relative eccentricities of the QC 
and the F orbits. The eccentricity of the QC run that we infer depends therefore on the 
residual eccentricity of the F run. A more intrinsic approach, used recently by Buonanno 
et al. [27], is to fit some eccentricity-dependent quantity to a full cycle (or more) of the 
orbital data. This approach yields similar, but somewhat smaller, eccentricities than 
those found here (despite our use of a QC orbit having larger initial separation and so 
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presumably smaller initial eccentricity). 
5.3. Waveform extraction 

We now turn our attention to the problem of extracting the gravitational wave signals 
from our numerical simulations using the Newman- Penrose quantity ^4. Given a spatial 
hypersurface with timelike unit normal n^, and given a spatial unit vector r M in the 
direction of wave propagation, the standard definition of ^4 is the following component 
of the Weyl curvature tensor, 

^4 = -C^^tm^m^ (27) 

where £^ = -^(n^ — r M ), and m M is a complex null vector (satisfying m M m M = 1) that is 
orthogonal to r M and n M . Here an overbar denotes complex conjugation. 

For (perturbations of) flat spacetime, \I/ 4 is typically evaluated on coordinate 
spheres, and in this case the usual choices for r M and are 

r» = (I)", (28a) 



^2r \d0 sin 9 d<p 

where (r, 9, <fi) denote the standard spherical coordinates. With this choice, ^4 can be 
expanded in terms of spin-weighted spherical harmonics of weight -2: 

4 (t, r, 9,<P) = J2 *T(t, r) - 2 Y lm (9, 0), (29) 

lm 



where the are expansion coefficients defined by this equation. 

For curved spacetime, there is considerable freedom in the choice of the vectors r M 
and m M , and different researchers have made different choices [27, 36, 37, 38, 39, 40, 7] 
that are all equivalent in the r — > 00 limit. We choose these vectors by first picking 
an extraction two-surface £ that is a coordinate sphere (r 2 = x 2 + y 2 + z 2 ) centered 
on the center of mass of the binary system (using the global asymptotically Cartesian 
coordinates employed in our code). We choose r^ to be the outward-pointing spatial unit 
normal to 8 (that is, we choose r, proportional to V«r) . Then we choose m M according to 
Eq. (28 ft), using the standard spherical coordinates 9 and <fi defined on these coordinate 
spheres. Finally we use Eqs. (27) and (29) to define the ¥^ coefficients. Note that 
our is not exactly null nor exactly of unit magnitude at finite r, so our definition of 
\I>4 m will disagree with the waveforms observed at infinity (and with those computed by 
other groups). Our definition does, however, agree with the standard definition given 
in Eqs (27)-(29) as r — > 00, so our definition only disagrees with the standard one by a 
factor of order 1 + 0(l/r). In this paper we compute ^4™ in the same way and at the 
same extraction radius for all runs, so the 0(1/ r) effects should not significantly affect 
our comparisons of these waveforms. 

Since our simulations use high spatial resolution all the way to the outer boundary, 
the outgoing radiation is fully resolved everywhere. Therefore, we could extract 
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Figure 7. Waveforms for the F2 run. Plotted are the six dominant coefficients, 
scaled by the factor IOOOtMadm- Solid lines represent the real parts and dashed 
lines the imaginary parts of The time axes are labeled in geometric units at the 
bottom, and in SI units for a 20+20 Mq binary at the top. 



waveforms at very large radii. The extracted wave signal lags the dynamics of the 
binary by the light-travel time to the extraction radius, and our evolutions currently 
fail shortly before merger. So extracting the wave signal at a very large radius would 
miss the most interesting part of the waveform close to merger. In order to retain most 
of the signal, we compromise by extracting the radiation at an intermediate distance: 
i?/Af adm — 57. Figure 7 presents the dominant waveform coefficients ^4. The 
coefficient is about a factor of ten smaller than the largest coefficient, ^f 2 . The \1/| 2 and 
\&4 6 coefficients are smaller by about another order of magnitude; and the \l/| 2 and ^f 4 
coefficients have amplitudes that are only about ~ 1/1000 that of ^if. 

5.4- Waveform comparisons 

In this section we make a number of quantitative comparisons between the waveforms 
produced by the evolution of quasi-circular, QC, initial data and those produced by the 
lower eccentricity, F, initial data. 

We can define a gravitational wave frequency associated with \I/ 4 m by writing 

*l m = A, m (*)e-*"»W (30) 

where Ai m (t) is its (real) amplitude and <f>i m (t) its (real) phase. The frequency, Qi m , 
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associated with is then defined as 

n lm = ^jp. (si) 

Figure 8 shows comparisons of the frequency of the dominant mode, ^22, from the QC 
and the F runs. This figure confirms the basic picture that emerged from our discussion 
in Sees. 5.1 and 5.2: a time offset At must be used to compare the QC and F runs 
properly; the QC run has an orbital eccentricity which causes Q22 to oscillate; and 
these oscillations are largely absent from the F run. Indeed, apart from the factor of 
two difference between orbital and the gravitational wave frequencies, the top panel of 
Fig. 8 looks very much like Fig. 2. This indicates that our coordinates are very well 
behaved — a feature that has also been observed in other numerical simulations, e.g. 
Ref. [41]. 

In order to make more detailed comparisons between the QC and the F waveforms, 
a phase offset A<p in addition to the time offset AT must be taken into account. These 
offsets are used then to redefine the waveform of the QC run: 

*i m Qc(*) = *4 m Q c(* - AT). (32) 

The same time and phase offsets are used for all values of I and m. Note that A0 
and AT represent differences between the QC and F evolutions. These offsets differ 
therefore from those often used in LIGO data analysis, where offsets are used to set the 
time and orbital phase at which a binary signal enters the LIGO band at 40Hz. 

We now estimate the phase offset Acf> needed in Eq. (32) to allow us to make direct 
comparisons between the QC and the F2 waveforms. We consider two effects: First, the 
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orbital phase of the time-shifted QC run differs from that of the F2 run by the phase 
accumulated by the F2 run during the time < t < AT. Second, the orbital frequencies 
of the QC and F2 runs differ, and this difference oscillates in time (cf. the right panel 
of Fig. 5), so the orbital phase difference between the two runs also oscillates in time. 
We take both of these effects into account: first, we evaluate the time-dependent phase 
difference, A0(t), between the waveforms of the time offset QC run, \&4Qc(t — At), 
and the F run, ^ 4F (t); second, we evaluate the time average of this A0(t) to obtain 
A(p ks 1.83. Using this value of Ac/> leads to waveforms for the QC and F2 evolutions 
that agree as well as can be expected in the presence of the other systematic errors, 
described below. 

The two gravitational wave polarizations, h+(t) and h x (t), are the real functions 
related to by 

\& 4 = h + -ih x . (33) 

Consequently, the - 2 Yi m components of h + (t) and h x (t) can be obtained by the double 
time integral, 

h ^(t) - ih l ™(t) = f dr [ T dr'¥ 4 m (r') + C lm + D lm t. (34) 

The constants Ci m and D lm account for the (unknown) values of h and h at the initial 
time tj. If the full waveform were known, they could be determined either at very early 
times or at very late times (i.e. after the merger and ringdown). Since we do not have 
complete waveforms for the present evolutions, we choose C\ m and Di m that make the 
average and the first moment of h l ™ x (t) vanish: 

2 dr h% (r) = = I' drr h% (r) . (35) 

The integration interval [ti,t 2 ] = [160Madm, 706Madm] is chosen to be the largest 
interval (excluding the initial transient radiation burst) on which data is available for 
both runs. 

Figure 9 shows the waveforms h l ™ for the evolution F2 (solid lines) and QC (dashed 
lines). To the eye, the waveforms look essentially identical. To quantify how well the 
two waveforms match, we use simple overlap integrals in the time domain: 

(36) 

where (hi,h 2 ) = dt h 1 (t)h 2 (t), and \\h\\ 2 = (h,h). The quantity /i gives the loss of 
signal to noise ratio obtained by filtering waveform hi with waveform h 2 - We evaluate 
the overlap integral in the time domain, rather than the frequency domain, to allow us 
to truncate the waveforms easily to the interval [ti,t2] during which both waveforms 
are available. During the evolutions presented here the gravitational-wave frequency 
changes by only a factor of two, so our decision not to weight by the LIGO noise 
spectrum should not change our results significantly for frequencies near the minimum 
of the noise curve. Furthermore, we evaluate \i directly for the different modes h l ™ x , 
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Figure 9. Waveforms h ™ (normalized by t/Madm) for the six dominant 
modes. The solid lines represent evolution of the low-eccentricity initial data (run F2). 
The dashed lines represent evolution of QC initial data time-shifted by AT = 59-Madm 
and phase-rotated by A<j> — 1.83. The time axes are labeled in geometric units at the 
bottom and in STunits for a 20+20 Mq inspiral at the top. 



Table 2. Waveform overlaps between the low-eccentricity run F2 and quasi-circular 
run QC (computed from runs with medium and high resolution). Each mode of QC 
has been time shifted and rotated by AT = 59A/adm and A<fi — 1.83. These numbers 
are subject to additional systematic effects as discussed in the text. 



mode 



nlh lm h lm 



high resolution 
) 



Ira i*lm ^ 
Fx: n QCxi 



medium resolution 



n QC+) 



L Fx ! "QCx, 



1=2, m=2 
1=3, m=2 
1=4, m=2 
1=4, m=4 
1=5, m=4 
1=6, m=6 



0.998 
0.997 
0.996 
0.991 
0.987 
0.981 



0.998 
0.997 
0.997 
0.991 
0.979 
0.980 



0.998 
0.997 
0.996 
0.993 
0.983 
0.986 



0.998 
0.998 
0.998 
0.993 
0.982 
0.982 



rather than for specific observation directions. This allows us to compare differences in 
the higher order modes with smaller amplitudes, which would otherwise be swamped 
by the dominant I — m — 2 mode. 

The overlaps between the QC and the F2 waveforms, obtained at AT = 59Madm 
and A<p = 1.83, are summarized in Table 2. Both medium and high resolution 
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overlaps are given in Table 2, confirming that the overlaps are not dominated by 
numerical errors. We note, however, that the medium resolution runs have more noise 
in the higher order modes at early times; so we shortened the integration interval to 
[ti, ^2] = [200Madm, 706Madm] to avoid contamination in those waveforms. 

The dominant uncertainty in the computed overlap /i arises because of our 
uncertainties in the integration constants Ci m and D\ m in Eqs. (34) and (35). Because the 
waveform has finite length, these constants are known only to an accuracy of ~ 1/N cyc , 
where A^ cyc is the number of cycles in the waveform. This error depends only on the 
length of the evolution, and can only be reduced by longer evolutions, not by higher 
resolution evolutions. We show in Appendix B (to lowest order in the uncertainties of 
Cim and Di m ) that the overlaps quoted in Table 2 are upper bounds. We also derive 
lower bounds for these overlaps there, which are smaller than the values given in Table 2 
by about 12/(7riV cyc ) 2 . So these lower bounds are about 0.02 smaller than the Table 2 
values for the m = 2 modes, and 0.002 smaller for the m = 6 modes. This systematic 
uncertainty is much larger than the mismatch of the waveforms for the m = 2 modes, 
so maximizing the overlaps by varying AT and A<p as independent parameters is not 
justified. 

6. Discussion 

In this paper, we have extended the quasi-equilibrium initial-data formalism to binary 
black holes with nonzero radial velocities. We have also used this formalism to construct 
initial data whose evolution results in very low eccentricity orbits: about an order of 
magnitude smaller than the orbits of quasi-circular initial data. 

The main differences between evolutions of the quasi circular, QC, and the low 
eccentricity, F, initial data sets are overall time and phase shifts: the QC initial data 
represents the binary at a point closer to merger. When we correct for these shifts, 
the orbital trajectories of the black holes and the gravitational waveforms they produce 
agree very well between the two runs. Various parameters measured in the QC run (e.g. 
orbital frequency or proper separation) oscillate around the corresponding values from 
the F run. The gravitational wave phase oscillates as well, but no significant coherent 
phase difference builds up during the five orbits studied here. We find waveform overlaps 
between the high-eccentricity and low-eccentricity runs of about 0.99. Therefore it 
appears that for the last five orbits before merger the differences between quasi-circular 
and low-eccentricity initial data are not important for event detection in gravitational 
wave detectors. Longer evolutions (e.g. equal mass binaries starting at larger separation, 
as well as unequal mass binaries with a longer radiation reaction time scale) have 
more cycles during which phase shifts could in principle accumulate. However, orbital 
eccentricity tends to decay during an inspiral and the orbital eccentricity in quasi- 
circular data should decrease as the initial separation increases, so longer evolutions are 
probably less sensitive to the eccentricity in the initial data. Thus we anticipate that 
the eccentricity of quasi-circular initial data will not play a significant role when longer 
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evolutions are used for event detection, but further study would be needed to confirm 
this. 

Finally, we note that construction of low-eccentricity inspiral initial data may be 
more difficult when the black holes carry generic spin. The process described in this 
paper merely adjusts the orbital parameters to obtain a trajectory without oscillations 
on the orbital timescale. For non-spinning equal-mass black holes sufficiently far from 
merger, a non-oscillatory inspiral trajectory seems to be a reasonable choice. But if 
non-negligible spins are present, this is not likely to be the case. For spins that are not 
aligned with the orbital angular momentum, the approximate helical Killing vector is 
lost, and there are likely to be a variety of oscillations on the orbital time scale. In 
these cases a more sophisticated model of the desired circularized orbit will be needed 
before a procedure for adjusting the orbital parameters to the appropriate values can 
be formulated. 
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Appendix A. Quasi-equilibrium initial data in inertial coordinates 

In this appendix we show that (ip C o, P l co ~ £rot, ct co ), where £* ot = (f2 x r)\ is a solution 
to the XCTS Eqs. (8a)-(8c) in asymptotically inertial coordinates (with appropriately 
modified boundary conditions) whenever (ip co , p i co ,a co ) is a solution in co-rotating 
coordinates. We also show that this solution leads to the same physical metric 7^ 
and extrinsic curvature as the original solution in co-rotating coordinates. The 
proof relies on three key observations: First, both solutions are assumed to make the 
same choice of free data Eqs. (9a) (96), (15a), and (156); second, the shift enters the 
XCTS equations and the boundary conditions (almost) solely through the conformal 
Killing operator, (L/3) J ; and third, £* ot is a conformal Killing vector, so (L£ rot ) u = 0. 
Hence the term —Q ot that is added to f3 l co (mostly) drops out of the equations. 

We first show that the XCTS equations remain satisfied: Since (L£ rot ) y = 0, it 
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follows from Eq. (5) that A i3 is unchanged by the addition of £* ot . So Eqs. (8a) and (86) 
remain satisfied. The only other shift containing term in Eq. (8 c) is f3 l diK, which 
vanishes because diK = from the choice of free data (K = 0) in Eq. (156); so Eq. (8c) 
also remains satisfied. 

We turn next to the boundary conditions: The boundary conditions used for the co- 
rotating coordinate representation of the XCTS equations are Eqs. (9c)-(9g) and (15c), 
while those used for the inertial frame representation are the same, except Eqs. (9 c?) 
and (9e) are replaced by Eqs. (16a) and (166). The boundary conditions, Eqs. (9c) and 
(15 c), depend only on ip and a and therefore remain satisfied. The apparent horizon 
boundary condition, Eq. (9/), implies the boundary condition on the conformal factor 
Eq. (13), which is unchanged since (L^ rot ) y = 0; and the new outer boundary condition, 
Eq. (16a), also holds because f3 l co satisfies Eq. (9d). 

The only remaining boundary conditions then are Eqs. (166) and {9g). Because 
= and (Tjj = 0, the null surface generated by k p coincides with the world tube of 
the apparent horizons, <Sah- The normal to this null surface is k^, because k p is normal 
to S by construction, and because k^k^ = 0. Therefore, in order for d t + CotA to be 
tangent to <Sah, as required by the boundary condition Eq. (166), it must be orthogonal 
to k p . The vector d t + C ot di has components an p + j3 p + £f ot , where (3^ = [0, (3 l \ and 
^ ot = [0,£ ot ]. Using k p = (n» + s»)/V2, it follows that 

= (d t + CoA) ■k=-^=[-a + (f3 i + . (A.l) 

This condition implies 

(3 i = as i -& ot + C on S, (A.2) 

with Qsi = 0, i.e., Eq. (17) in the main text. So the boundary condition Eq. (166) is 
satisfied because (3 l co = as 1 + ( l satisfies Eq. (12). 

The vector £ l that appears in Eq. (A.2) is further constrained by the shear boundary 
condition, Eq. (9g), which we consider next. The shear is defined as 

o~fj,u = -L AtJy p ' Vpfcff, (A. 3) 

where _L po = h^ p h u a ^ — \h tliV h po ' . Substituting Eq. (11) into this expression, and 
subsequently using Eqs. (2), (4), and (5) results in 

V2a tJ = --j- V [(L/3) w - ^u kl ] + V V fc s,. (A.4) 

For any vector field v l decomposed into normal and tangential parts, v l = v m s m s l 
it follows that 

± tJ kl (hv) kl = (hsv\\)ij + 2v m s m L i3 kl V kSl . (A.5) 
Using this identity and Eq. (17), the shear can be rewritten as 

V2a tJ = -L L i3 kl [(hUh + ^u kl ] - ^(LsCk- (A.6) 
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Once more, £* ot drops out because it is a conformal Killing vector. Also, since w^- = 
by Eq. (9 a), we find that the shear vanishes iff C is a conformal Killing vector within 
the 2-surface S: 

a lj = = (h s (f. (A.7) 

Equation (18) now follows from the identity (L s (y j = ^ _4 (L,s£) u . This implies then 
that the boundary condition Eq. (9g) is satisfied since it is assumed to be satisfied in 
the co-rotating case. 

Finally, we note that the physical metric 7^ and extrinsic curvature produced 
by the inertial frame version of the problem are identical to those of the original co- 
rotating frame version. The conformal metric and conformal factor are identical in the 
two versions, so the physical metrics are identical trivially from Eq. (3). Since £* ot is a 
conformal Killing vector, it follows that is identical from Eq. (5); so it follows from 
Eq. (4) (with K — 0) that the extrinsic curvatures are identical as well. 



Appendix B. Errors caused by finite-length waveforms 

The error in the waveform overlaps caused by the uncertainty in the integration 
constants can be determined as follows: Denote our numerically computed waveforms 
by h x + s x , where h x stands for the unknown "true" waveform obtained with the correct 
values of the integration constants, and e x represents the error introduced by computing 
these constants with a truncated waveform. The label x stands for either F or QC. 
The quantity of interest is the overlap between the "true" waveforms, 

Mft F ,ft 0C )= iif'' l y l i . (b.i) 

\\iIf\\ \\iiqc\\ 

where (h u h 2 ) = f£ h^fatydt, and \\h\\ 2 = (h,h). The errors e x are those caused 
by the uncertainty in the constants Ci m and Di m in Eq. (34), and the e x are therefore 
linear functions of time. Furthermore, choosing the integration constants by Eq. (35) 
makes the numerical waveforms h x + e x orthogonal to functions linear in time, so that 
(h x + e x , e y ) = 0, where x, y G {F, QC}. Using this result, and neglecting terms of order 
0(e 3 ), one finds 

n(h F + e F , h QC + e QC ) = fJ,(h F , h QC ) 

j. ,.(h h a ( W £f W" 1 W £ Qc\\ 2 ( £ F,e QC ) \ (n9) 

It is straightforward to show that /i(hf,hQc) = 1 — 0(5h 2 ) where 5h = h F — hqc- 
Therefore, replacing fj,(h F , hQc) — > 1 in the last term of Eq. (B.2) changes the result 
only by terms of order 0(5h 2 e 2 x ). Furthermore, replacing |]/iqc|| ~~ * W^fW m the 
denominators of Eq. (B.2) affects the result only by terms of order 0(5he 2 ). Neglecting 
both of these higher order contributions, we find 



nih F + e F , h QC + e Q c) = fJ>(h F , h QC ) + — oM , M2 — • (B.3) 



\\ £ F 


_ £ Qc\ 


12 


2| 


\h F \ 


2 
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\e F 


-£q\ 


I 2 / (IM 


+ lk<3cll) 2 


2| 


\h F \ 


2 


2 




2 



Because the last term is non-negative, the "true" overlap /i(h F , Hqc) is always smaller 
than the numerically computed overlap fj,(h F + £F,hqc + £qc)- Using the triangle 
inequality, we can bound the last term in Eq. (B.3) by the error | \e x \ | 2 /| \h x \ | 2 in either 
the F or the QC waveform: 

2^£- (B-4) 

\\h x \\ 

Finally, we estimate | \e x \ | 2 /| \h x \ | 2 by applying Eqs. (34) and (35) to a pure sine- 
wave: h(t) = sin(t). It is straightforward to evaluate the integrals in Eq. (35) for this 
simple case, giving the bound | |e| | 2 /| \h\ | 2 < 6/(7riV cyc ) 2 , where A^ cyc = (t 2 — t 1 )/(2n) is 
the number of cycles in the interval [ti,^]- Therefore, we arrive at the bounds 

12 

H{h F + e F , h QC + e Q c) > n(h F , h QC ) > n{h F + e F , h QC + e QC ) T*j2~> ^ B - 5 ) 

^ eye 

as mentioned in the main text. 
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